*** Set working directory ***
cd "C:\..."

*********************** Load Data **********************************************
clear
graph drop _all
cap drop all

import excel "Data.xlsx", sheet("data") firstrow

************************** Setup ***********************************************
* Choose impulse response horizon
local hmax = 36

* Time variable
gen monthly_date = mofd(date)
format monthly_date %tm
drop date
rename monthly_date date
sort date
tsset date
egen time = group(date)
drop if date == .

* Standardize MP Shock series - unit standard deviation
egen mps_std = sd(chg_sveny02)
gen mps = chg_sveny02/mps_std

sort date
tsset date

* Normalize GSCPI for 1998-2019 sample
rename gscpi gscpi_unscl
egen gscpi_mean = mean(gscpi_unscl)
gen gscpi_dm = gscpi_unscl - gscpi_mean
egen gscpi_std = sd(gscpi_dm)
gen gscpi = gscpi_dm/gscpi_std

* GSCPI mp interaction terms
gen gscpi_mps = L.gscpi*mps
gen gscpi_mps1 = L.gscpi*L.mps
gen gscpi_mps2 = L.gscpi*L2.mps

* Put variables in log form
gen cpi_inf = ln(cpi)
gen rgdp_growth = ln(ip)
rename unrt unrt_chg
rename retail retail2
gen retail = ln(retail2)

* SP500 in log form
gen stk = ln(sp)


* Endog variables
forvalues h = 0/`hmax' {
	gen onetp`h' = f`h'.onetp
	gen twotp`h' = f`h'.twotp
	gen fivetp`h' = f`h'.fivetp
	gen tentp`h' = f`h'.tentp
}

************************** LP regressions **************************************

* 1 year TP
eststo clear
cap drop b1 u1a d1a u1 d1 Years1 Zero1
gen Years1 = _n-1 if _n<=`hmax'
gen Zero1 =  0    if _n<=`hmax'
gen b1=0
gen u1a=0
gen d1a=0
gen u1=0
gen d1=0

forv h = 0/`hmax' {
	* levels
	newey onetp`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(0/2).rgdp_growth l(0/2).cpi_inf l(0/2).unrt_chg l(0/2).retail l(1/2).onetp, lag(`h')
replace b1 = _b[gscpi_mps]                    if _n == `h'+1
replace u1a = _b[gscpi_mps] + 1.645* _se[gscpi_mps]  if _n == `h'+1
replace d1a = _b[gscpi_mps] - 1.645* _se[gscpi_mps]  if _n == `h'+1
replace u1 = _b[gscpi_mps] + 1* _se[gscpi_mps]  if _n == `h'+1
replace d1 = _b[gscpi_mps] - 1* _se[gscpi_mps]  if _n == `h'+1

eststo
}


* 2 year TP
cap drop b3 u3a d3a u3 d3 Years2 Zero2
gen Years2 = _n-1 if _n<=`hmax'
gen Zero2 =  0    if _n<=`hmax'
gen b3=0
gen u3a=0
gen d3a=0
gen u3=0
gen d3=0

forv h = 0/`hmax' {
	* levels
	newey twotp`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(0/2).rgdp_growth l(0/2).cpi_inf l(0/2).unrt_chg  l(0/2).retail l(1/2).twotp, lag(`h')
replace b3 = _b[gscpi_mps]                    if _n == `h'+1
replace u3a = _b[gscpi_mps] + 1.645* _se[gscpi_mps]  if _n == `h'+1
replace d3a = _b[gscpi_mps] - 1.645* _se[gscpi_mps]  if _n == `h'+1
replace u3 = _b[gscpi_mps] + 1* _se[gscpi_mps]  if _n == `h'+1
replace d3 = _b[gscpi_mps] - 1* _se[gscpi_mps]  if _n == `h'+1

eststo
}


* 5 year TP 
cap drop b5 u5a d5a u5 d5 Years3 Zero3
gen Years3 = _n-1 if _n<=`hmax'
gen Zero3 =  0    if _n<=`hmax'
gen b5=0
gen u5a=0
gen d5a=0
gen u5=0
gen d5=0

forv h = 0/`hmax' {
	* levels
	newey fivetp`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(0/2).rgdp_growth l(0/2).cpi_inf l(0/2).unrt_chg  l(0/2).retail l(1/2).fivetp, lag(`h')
replace b5 = _b[gscpi_mps]                    if _n == `h'+1
replace u5a = _b[gscpi_mps] + 1.645* _se[gscpi_mps]  if _n == `h'+1
replace d5a = _b[gscpi_mps] - 1.645* _se[gscpi_mps]  if _n == `h'+1
replace u5 = _b[gscpi_mps] + 1* _se[gscpi_mps]  if _n == `h'+1
replace d5 = _b[gscpi_mps] - 1* _se[gscpi_mps]  if _n == `h'+1

eststo
}


* 10 year TP
cap drop b7 u7a d7a u7 d7 Years4 Zero4
gen Years4 = _n-1 if _n<=`hmax'
gen Zero4 =  0    if _n<=`hmax'
gen b7=0
gen u7a=0
gen d7a=0
gen u7=0
gen d7=0 

forv h = 0/`hmax' {
	* levels
	newey  tentp`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(0/2).rgdp_growth l(0/2).cpi_inf l(0/2).unrt_chg l(0/2).retail l(1/2).tentp, lag(`h')
replace b7 = _b[gscpi_mps]                    if _n == `h'+1
replace u7a = _b[gscpi_mps] + 1.645* _se[gscpi_mps]  if _n == `h'+1
replace d7a = _b[gscpi_mps] - 1.645* _se[gscpi_mps]  if _n == `h'+1
replace u7 = _b[gscpi_mps] + 1* _se[gscpi_mps]  if _n == `h'+1
replace d7 = _b[gscpi_mps] - 1* _se[gscpi_mps]  if _n == `h'+1

eststo
}


************************* Plot IRF's *******************************************
* 1 year TP
twoway ///
(rarea u1a d1a  Years1,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u1 d1  Years1,  ///
fcolor(gs12) lcolor(white) lpattern(solid))  ///
(line b1 Years1, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero1 Years1, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white)) ///
title("One-Year Term Premium", size(small)) ///

gr rename fig_1tp_int, replace


* 2 year TP
twoway ///
(rarea u3a d3a  Years2,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u3 d3  Years2,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b3 Years2, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero2 Years2, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white)) ///
title("Two-Year Term Premium", size(small)) ///

gr rename fig_2tp_int, replace


* 5 Year TP
twoway ///
(rarea u5a d5a  Years3,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u5 d5  Years3,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b5 Years3, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero3 Years3, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white)) ///
title("Five-Year Term Premium", size(small)) ///

gr rename fig_5tp_int, replace


* 10 year TP
twoway ///
(rarea u7a d7a  Years4,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u7 d7  Years4,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b7 Years4, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero4 Years4, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white)) ///
title("Ten-Year Term Premium", size(small)) ///

gr rename fig_10tp_int, replace


* Combine all 4 IRFs
gr combine fig_1tp_int fig_2tp_int fig_5tp_int fig_10tp_int, ///
rows(2) cols(2) graphregion(color(white)) plotregion(color(white)) ///

********************************************************************************